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The dynamical behavior of a harmonic chain in a spatially 
periodic potential (Frenkel-Kontorova model, discrete sine- 
Gordon equation) under the influence of an external force 
and a velocity proportional damping is investigated. We do 
this at zero temperature for long chains in a regime where 
inertia and damping as well as the nearest-neighbor interac- 
tion and the potential are of the same order. There are two 
types of regular sliding states: Uniform sliding states, which 
are periodic solutions where all particles perform the same 
motion shifted in time, and nonuniform sliding states, which 
are quasi-periodic solutions where the system forms patterns 
of domains of different uniform sliding states. We discuss the 
properties of this kind of pattern formation and derive equa- 
tions of motion for the slowly varying average particle density 
and velocity. To observe these dynamical domains we suggest 
experiments with a discrete ring of at least fifty Josephson 
junctions. 

PACS numbers: 74.50.+r, 46.10.+Z, 03.20.+i 



I. INTRODUCTION 

60 years ago, Frenkel and Kontorova introduced a sim- 
ple model which has become popular in many fields of 
solid-state physics and nonlinear dynamics jlj. They in- 
vented their model in order to describe the motion of 
a dislocation in a crystal ||. Meanwhile, the Frenkel- 
Kontorova (FK) model has become also a model for an 
adsorbate layer on the surface of a crystal for ionic 
conductors ||D, f° r glassy materials |5],[4].[5j! for charge- 
density- wave (CDW) transport ||], for chains of coupled 
Josephson junctions ^,^|, and for sliding friction pj. 

The FK model is a chain of particles with mass m 
coupled by a harmonic nearest-neighbor interaction with 
stiffness n. ft is under the influence of an external spa- 
tially periodic potential with periodicity c and strength 
Uq. Here we study the nonequilibrium behavior of the 
FK model driven by a force F. We assume energy dis- 
sipation due to a usual damping force with a damping 
constant 77. After rescaling time and space, one gets the 
following equation of motion in dimcnsionless units: 



implies that the average particle density l/a is constant, 



i.e. 



1 — 2xj — bshiXj 



F. 



(1) 



b = (2tt/c) 2 U /k, and F = 



where 7 = n/\ 
(2it/c)F/k. The time and space units are \Jmj k and 
c/27r, respectively. We assume periodic boundary condi- 
tions, i.e., 
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a = 2tt ■ 
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where N is the number of particles and M is an arbi- 
trary integer. Note that the periodic boundary condition 



Due to symmetries a can be restricted to a S [0, n] with- 
out loss of generality. 

The dynamical behavior of the FK model has already 
been studied in several limits: (i) in the overdampcd limit 
(i.e., ij can be neglected) for large N [| 11 1, (ii) in the 
limit of zero damping and driving (i.e., 7 = F — 0) for 
a/2-K near integer values allowing well separated 2tt kinks 
|12[ as well as for the most incommensurate value (i.e., 
the golden mean) of a/27r JDjfl , and (iii)in the under- 
damped case for small N (N < 10) |,fij. In a series 
of papers Braun et al. have studied recently the un- 
derdamped dynamics of a generalized FK model with 
N > 100 but a near zero and ir p^!7[|. 

In this series of two papers, we study the underdamped 
FK model for large numbers of particles (i.e., N > 100). 
We do not restrict ourselves to values of a/2ir near inte- 
ger or half integer values where the dynamical behavior 
can be described in terms of kinks. The system has sta- 
tionary, periodic, quasi-periodic, and chaotic solutions. 
Of special interest is the transition from stationarity to 
sliding, the so-called pinning- depinning transition, and 
the backward process. In the first paper, we investi- 
gate the periodic and quasi-periodic solutions. The sec- 
ond paper will be concerned with the depinning-pinning 
transition between stationary states and spatio-temporal 
chaos. Preliminary results have already been published 
in a conference proceedings . 

In section IV we will see that the FK model sponta- 
neously forms spatial-temporal patterns as many other 
spatially extended systems driven far from equilibrium 
p9| . These patterns are caused by bistability and in- 
stability of the uniform sliding state. In the uniform 
sliding state all particles perform exactly the same reg- 
ular and periodic motion. Different particles differ only 
in the phase of this motion. The phase difference of two 
neighboring particles is everywhere the same. That is, 
the density of particles is on average constant along the 
chain. The uniform sliding state is the only nonstation- 
ary state in the overdamped limit f2(i|| but it has also be 
studied in the underdamped case |9|j21[|. 

The discreteness of the chain leads to several reso- 
nances in the underdamped limit ||Jl^,^l|,^|. The con- 
sequence is bistability. For spatially extended bistable 
systems it is well-known that domain-like patterns are 
possible [fl9|| . In the FK model these domains can be 
characterized by the average particle density 1 /a and the 
average particle velocity v. We find that states with two 
or three different types of domains survive in the long- 
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time limit. The number and the width of the domains 
can vary, leading to a quite large number of possible so- 
lutions. Because of different average velocities in differ- 
ent domain types, the motion in the domain solutions is 
quasi-periodic. Assuming that a and v are slowly vary- 
ing functions in time and space, we derive an approx- 
imate equation of motion for them. With its help we 
understand why the domains do not disappear in the 
long-time limit and why there are not more than three 
different domain types. Furthermore it turns out that 
the state in the long-time limit can be understood as a 
spatially chaotic solution of a corresponding dynamical 
system. A special variant of these domain solutions is 
the traffic-jam state where the particle velocity in one 
domain is zero |T^-|l7|]. 

The paper is organized as follows: In Section |d] we de- 
rive two different but mathematically equivalent formulas 
for the relation between the force F and the average slid- 
ing velocity v. In Section III the periodic solution (i.e., 
the uniform sliding state) and its stability are discussed. 
The domain- like states are investigated in Section IV. 
In Section ^ the main part of the paper concludes with 
some remarks concerning possible experimental observa- 
tions of these domain-like states and similarities to other 
pattern forming systems. The Appendices describe our 
numerical and analytical scheme to obtain the uniform 
sliding state and to analyze its stability. 



II. AVERAGE SLIDING VELOCITY AND 
EFFECTIVE FRICTION FORCE 

The average sliding velocity v of the chain reads 

r J , N 

m — / 

T- 



V 



(4) 



Plotting the measured or calculated values of v for dif- 
ferent values of the applied force F one gets the so- 
called velocity-force characteristic. In CDW systems and 
Josephson-junction arrays it corresponds to the current- 
voltage and voltage-current characteristics, respectively. 
In the context of friction, F can be interpreted as the ef- 
fective kinetic friction as a function of the sliding velocity 
v. 

There are two mathematically equivalent relationships 
between F and v. The first one can be obtained by taking 
the time average of the sum of ([!]) over all particles: 
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8XD.X j(t) dt. 



(5) 



This formula assumes that the average acceleration of the 
chain is zero. The other relationship can be derived from 
the fact that the energy released during sliding has to be 
dissipated totally, i.e., 
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(6) 



This result clearly shows that in the case of no dissipation 
(i.e., 7 = 0) the chain may slide without any applied 
force. In Section [III_B| we will discuss the condition of 
this possibility. Eq. (|6|) also shows that the mobility v / F 
is always less than or equal to I/7, i.e., the mobility in 
the limit 6^0. 



III. UNIFORM SLIDING 

Because of the symmetries of the equation of motion 
there exist nonstationary solutions called uniform sliding 
states. They are characterized by the fact that every 
particle performs the same motion but shifted in time. 
That is, x J+ i(t) = Xj(t + Ti), for j = 1, . . . , N. Thus we 
need only one function, the dynamic hull function /, to 
describe the motion of all particles |2l| Jll|,|8|,p[ : 



Xj(t) = ip + aj + vt + f(tp + aj + vt), 



(7) 



where v is the average sliding velocity and ip is an arbi- 
trary phase. Because of the discrete translation symme- 
try of (0) the hull function is periodic, i.e., 



/(p + 27r) = /(p). 



(8) 



Plugging the ansatz (Q) into the equation motion (|l|) 
loads to a differential delay equation for the hull function 

/(¥>): 

v 2 f" (ip) + jv[l + f'(<p)] = f(ip + a) + f(<p -a)- 2/fo>) 

-bw[<p + f(<p)]+F. (9) 

The average sliding velocity v has to be chosen in such 
a way that a 27r-periodic function fulfills (||). If /(<£>) 
is a solution then f(<p + ip)+ip is also a solution. Wo 
make the definition of the dynamic hull function unique 
by restricting the solutions of (^) to 27r-periodic functions 
with zero average, i.e., 



2tt 



f(<p) d V = 0. 



(10) 



For the uniform sliding state @, the relationships 
and (H) between v and F reduce to 



and 
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i[tp + f(<p)]d(p 



(11) 



(12) 



Using (Bl) we get: 



respectively. 

Instead of solving (||) for agiven value of F it is more 
convenient to replace F by ([□]) or ([l2|) and solve @ for 
a given value of v. From the solution f(<p) one obtains 
the corresponding F. In the overdamped limit, it is well- 
known that F is a monotonically increasing function of 
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v 1 20 §]. In the underdamped case resonances lead to 
nonmonotonic velocity- force characteristics J23|. 
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FIG. 1. The dynamic hull function f(ip). The inset shows 
the corresponding particle position as a function of time. The 
parameters are a/2iv = (3 — yE)/2, b — 2, v — 0.75, and 
7 = 0.5. Because the solution is obtained by the method 
described in Appendix |a] the average sliding velocity v is pre- 
scribed. The corresponding value of F is approximately 0.74. 
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FIG. 2. Velocity-force characteristic of the uniform slid- 
ing state ((?]). Solid (dotted) lines indicate stable (unstable) 
solutions. Dashed and dashed-dotted lines in the inset d enot e 
analytically obtained approximations given by ( |l5| ) and (A9), 
respectively. The parameters are a/2-K = (3— v5)/2, 7 = 0.5, 
and 6=1 (left curve) and b = 2 (right curve and inset). The 
arrows denote the resonant values of v given by @. The 
numbers indicate the order n. 



A. Resonances 

We solve the hull-function equation (^|) numerically by 
expanding / into a Fourier series. The details are de- 
scribed in Appendix |a|. Figure [l] shows an example of a 
dynamic hull function. 

In the underdamped case the numerically obtained 
velocity-force characteristics exhibits clearly peaks as can 
be seen in figure ||. Near these peaks an increase of the 
driving force F by a considerable amount leads only to a 
slight increase of the average velocity v. In other words, 
the differential mobility, dv/dF, is much lower than the 
mobility without any periodic potential (i.e., I/7). The 
reason for this behavior is that the additional energy is 
only partially turned into a larger kinetic energy of the 
center of mass of the chain, whereas the main part is 
turned into oscillatory motion of the particle due to res- 
onances. 

In the frame of the center of mass which moves with av- 
erage sliding velocity v, the external potential leads to a 
time-periodic force acting on the particle. The frequency 
of this force, the so-called "washboard frequency", is 
given by the velocity of the center of mass divided by the 
period of the potential. Resonance occurs if the wash- 
board frequency is equal to the eigenfrequency of the 
phonon with wave number k = a. To see this, we solve 
(0) in the approximation sin[tp + f((f)] ~ simp. We get 



/O) 



ib/2 



v + ijv 



where to(k) is the phonon dispersion relation: 



Lu(k) = 2 



sm 



(13) 



(14) 



In order to get the corresponding value of F one can 
use either (11) or (|l2|). Although both equations are 
equivalent, we get for the approximation ( |l3| ) different 
results. Evaluating ([fl]) leads to the obviously wrong 
result F — jv. This is can be understood from the fact 
that ( p"3| ) is only the leading term of an expansion of 
the hull function in powers of b. Thus F — 71; + 0(b 2 ). 
Instead of calculating the next order in / one can use 
(O) which leads to 



F = jv 1 + 



b 2 /2 



[ W *(o) 



.,212 



+ 7^u 



2/i.2 



0(b 4 ) 



(15) 



Note that F has to be an even function of b because 
the external potential b cos x is an odd function of a; — 
7r/2. For large values of v, ( |l3| ) approaches zero and 
therefore F — > jv. That is, if the washboard slides very 
fast, the particles cannot follow the fast pushing by the 
washboard. Thus the chain slides like a rigid solid. We 
call this state the solid-sliding state. 

Figure || shows that for large values of b (or small val- 
ues of 7) this simple approximation (|l^) overestimates 
the strength of the resonance line. That is, there is an 
effective damping, larger than 7, which increases with in- 
creasing oscillation amplitude. This larger damping can 
be understood by phonon coupling due to the nonlin- 
earity in the equation of motion. This coupling opens 
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up additional channels for energy dissipation yielding a 
higher effective damping constant. The inset of Fig. ^ 
shows that a remarkably good approximation of these 
additional dissipation processes is given by the Galerkin 
approximation f(ip) = A cos(ip + ip). Projecting the hull- 
function equation (|9]) together with this ansatz onto sin (p 
and cos ip leads to two transcendental equations for A 
and i\>. After elimination of i\> one can parameterize 
the velocity-force characteristic by the amplitude A [see 
Eq. ([A|) in Appendix |A|]. 

In order to understand the other resonance peaks seen 
in Fig. |]one has to keep in mind that the external poten- 
tial not only leads to a periodic driving force but also to 
a modulation of the eigenfrequencies. Thus parametric 
processes make it possible to excite other phonons. In the 
framework of a perturbation theory with b as the small- 
ness parameter, the elementary processes corresponding 
to these resonance lines are the decay of n "washboard 
waves" (wavenumber a, frequency v) into a single phonon 
with wavenumber k and frequency u(k). Assuming mo- 
mentum and energy conservation one gets k = na and 
v = v n , where v n is the superharmonic resonance fre- 
quency of order n: 



uj(na) 



(16) 



The positions of the first few resonances are shown in 
Fig. ||. The agreement with the actual positions of the 
resonance peaks is quite good. Near the superharmonic 
resonance of order n the n-th Fourier mode of the hull 
function has a maximum. For example, the square in 
Fig. ^ corresponds to the hull function in Fig. [j] which 
is clearly dominated by exp(2itp). Because of finite dis- 
sipation, superharmonic resonances of higher order may 
be hidden behind a nearby resonance of lesser order as, 
e.g., the fourth resonance in Fig. pi 

The superharmonic resonance has been already inves- 
tigated in the literature, experimentally as well as 
theoretical ly Q ,[il,^l|,p2| . The superharmonic resonance 
condition (flq) was observed first in numerical experi- 
ments by Aubry and de Seze |^l|,^|. They studied the 
FK model without damping but with a very small driving 
force. They found that the velocity of the center of mass 
did not increase linearly in time but it was locked for fi- 
nite time intervals at velocities given by ([l6]). They also 
studied at the first time the underdamped and driven FK 
model J2l| . Peyrard and Kruskal found the same locking 
phenomenon for the velocity of a 27r-kink fl2|| . In this 
case the resonance frequencies are given by 

- nS w(no) with cj{k) = Jb + ^Hk). (17) 
n 

Ustinov et al. also observed resonance lines in numeri- 
cal simulations of the damped and driven FK model [jl4| . 
Van der Zant et al. have reported evidence of these res- 
onances in experiments with a ring of eight Josephson 
junctions p4j . 

Ustinov et al. as well as van der Zant et al. explained 
the found resonances by the following mechanism which 
leads to (|l7|). The mechanism relies on the assumption 
that 27r-kinks are traveling in the ring. Most of the time 



the particles are in a potential well. When a kink travels 
through a particle it jumps into the neighboring well and 
oscillates. Because of the periodic boundary conditions 
the jumps occur in equidistant time steps. The distance 
between two kinks in terms of the number of particles in- 
between is given by N/M — 2tt /a. The kink velocity c is 
related to the average sliding velocity v by v — ac. Super- 
harmonic resonances of order n occur if the time interval 
between two jumps, i.e., (N/M)/c = 2tt/v is n times the 
oscillation time of the particles, i.e., 2nn/uj{k) 1 where 
u>(k) given by (|l7j ) is the dispersion relation of the lin- 
earized equation of motion (Q). The wavenumber k times 
the distance N/M has to be 2irn, that is, k — an. Ac- 
cording to this mechanism the superharmonic resonance 
condition is therefore (|lj). This picture is valid only if 
the distance between two kinks, i.e., N/M — 2ir/a, is 
much larger than one. That is, the motion can be de- 
scribed by kinks only if a/27r is near an integer value. 
Thus, we expect that (0) is valid for a = 0(1) and fll~7j ) 
for a <C 1. For example, the numerically obtained values 
of positions of the resonance peaks reported by Watan- 
abe et al. || are more close to (|l^ ) than to ( fl7P because 
a > 1. 



B. Instabilities 

In order to discuss the instability of the uniform sliding 
state (Q), one has to investigate the dynamics of small 
perturbations Sxj . They are governed by the equation of 
motion (P linearized around (Q): 

5xj + jSij — Sxj-i + Sxj+i — 2Sxj 

—bcos[aj + vt + f(aj + vt)]Sxj. (18) 

The periodic boundary condition (|^) turns into 5xj+ n — 
Sxj. In accordance with the Floquet-Bloch theorem one 
can write any solution of this equation as a sum of solu- 
tions of the form 



5xj(t) — Ck(aj + vt)e N 



1, 



,N, (19) 



where Ck{f) is a 27r-periodic function and is the so- 
called Floquet exponent. The uniform sliding state is 
stable if the real part of Afc is negative for all values 
of k. There is always a solution with Ao = 0. It is 
the Goldstone mode 5xj = dtXj = v + vf'{aj + vt) = 
co(aj + vt), which follows directly from differentiating 
(^). Appendix [b| describes the scheme we have used to 
solve (18|) with the Floquet-Bloch ansatz (|l9| ) numeri- 
cally. The dotted lines in Fig. ^ denote unstable parts of 
the velocity-force characteristics. 

Two different mechanisms may lead to an instability 
of the uniform sliding state. The first one is negative 
differential mobility, i.e., a negative slope in the velocity- 
force characteristic. A small positive velocity fluctuation 
v — > v + Sv accelerates the chain because the applied 
force is larger than the force necessary to keep the new 
velocity v + Sv constant. 

The second type of instabilities is caused by parametric 
resonance. In the framework of multi-phonon process 
it corresponds to the decay of n washboard waves into 
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two phonons with wavenumber na/2 ± q. Parametric 
resonance of order n can be expected for values of the 
average sliding velocity v which are given by 



Wfront 



Vn(9) = 



Lo(na/2 + q) + uj(na/2 — q) 



(20) 



Because parametric resonance is a threshold phe- 
nomenon, the amplitude b of the washboard wave has 
to exceed a critical value which is proportional to 7 1 /™ 
fl9| , p5f . This is true only for values of v between the min- 
imum and the maximum of (pp|). For velocities outside 
this interval parametric resonance is still possible but the 
threshold increases. For zero damping the uniform slid- 
ing state is unstable for any value of v below a certain 
critical value v c which is approximately calculated in Ap- 
pendix |b| 



v c 



16 cos 2 



2b. 



(21) 



The actual value of v c obtained from the numerical sta- 
bility analysis agrees very well with this formula even for 
large values of b. The numerical value of v c is less than 
( pl| ) but deviates not more than 10 percent for b < 4. 
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FIG. 3. Dynamical domains of different particle densi- 
ties. A series of snapshots taken at equidistant time steps 
(St = 2n/v) are shown. Each dot denotes the position of a 
particle. A particular zig-zag shaped snapshot is highlighted. 
The zigs and zags correspond to two different kind of domains 
that are characterized by uniform sliding. The parameters are 
N = 144, M = 55, b = 2, 7 = 0.5, and F = 0.8. 



IV. SLIDING DOMAINS 

Near resonance peaks the system is bistable and has 
therefore the opportunity to organizes itself into domains 
of different uniform sliding states. Figure || shows a typi- 
cal example with ten domains. There are only two types 
of domains. Each domain is characterized by uniform 
sliding. That is, in each domain the particle motion is 
given by ([?]) but the hull function / and the value of a 
and v are different in each domain. One can say that 
the domains are characterized by different particle densi- 
ties 1/a. Neighboring domains are separated by domain 
walls (fronts) of finite size. Conservation of the number 
of particles implies that a front has to move with the 
velocity 



a 2 vi - aiv 2 
a 2 — a\ 



(22) 



where the average particle distance and the average par- 
ticle velocity of each domain type is given by 01,2 and 
i>i,2) respectively. From the viewpoint of the particles, 
we can express the front velocity in terms of how fast the 
front travels from one particle to the next. It is given by 



a 2 - cii 



(23) 



Because (|22|) and (|23|) are symmetric in the exchange of 
the indices, all fronts propagate with the same velocity 
leaving the widths of the domains constant. The num- 
bers A^i^ of particles in each type of domain fulfill the 
following constraints: 

Ni+ N 2 = N and aiiVi + a 2 N 2 = 2wM = aN (24) 

Because of < Ni, 2 < N the particle density for one 
type of the domains is larger than 1/a whereas for the 
other type it is less than 1/a. In the following the type 
with the larger density will be number one. Thus, 



ai < a < a 2 . 



(25) 



The average sliding velocity v of a domain state is given 
by 



V1N1 + v 2 N 2 
N ' 



(26) 



A domain-type state is in general a quasiperiodic mo- 
tion with three frequencies: v\ and v 2 from the periodic 
motion in each domain type and 2irc/N from the cyclic 
motion of the fronts through the system. 
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FIG. 4. (a) The velocity-distance characteristic for the 
uniform sliding states. Solid (dotted) lines indicate stable 
(unstable) states. The dash-dotted horizontal (tilted) line de- 
notes v — F/ry (the main resonance vi). The open (filled) 
squares denote the initial (final) domain states of the simula- 
tion shown in Fig. ^. (b) The numerically obtained values of 
the average particle distances aj behind the front for a given 
value at in front of the front. Those data points are connected 
by a solid line where a continuous function is expected. The 
dotted line denotes the inverse function. The circle near the 
intersection of both functions denote the values of ai and 02 of 
the numerically found two-domain solution. The parameters 
are b = 0.5, 7 = 0.5, and F = 0.6. 
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A. The two-domain state 

To understand this kind of pattern formation we inves- 
tigate in detail two-domain states. First, we have to dis- 
cuss the relationship between the average sliding velocity 
v and the average particle distance a for the uniform slid- 
ing state at a fixed value of F, i.e., the velocity- distance 
characteristic. A typical example is shown in figure ^(a). 
As in the velocity-force characteristic, resonances are re- 
sponsible for folds. 
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FIG. 5. Snapshots of the evolution of a two-domain state 
taken at equidistant time steps (At = 12). The presentation 
is the same as in Fig. |j| The initial state is a two-domain state 
with ai = 27r/100 and a? = 2n/5. The initial positions and 
velocities of the particles are calculated by using (Q). The 
parameters are N = 250, M = 25, b = 0.5, 7 = 0.5, and 
F = 0.6. 

A two-domain state is completely characterized by two 
points on the stable branches of the velocity-distance 
characteristic which fulfill (^5|). The velocity c of the 
front is the slope of the line connecting both points (see 
Fig. [7]) . The sizes of the domains are determined by the 
solutions of (^) . 

From this consideration one would expect a continu- 
ous family of two-domain states parameterized by two 
real numbers. But this is not the case, as figure ^| clearly 
shows. In our numerical simulations we always found 
that the system selects dynamically the same pair of 
points on the velocity- distance characteristic. This is 
true even if a is changed as long as a £ (ai,a 2 ). How 
does the system selects a certain pair of values for a^? 
A careful inspection of Fig. || reveals that behind the 
fronts new domain states are selected. These states are 
independent of the initial states. The interface between 
the new domain state and the old one does not form 
a front. It smears out and in the long-time limit the 
domain states approach to uniform densities with well- 
defined values of a. Numerical experiments show that 
the state af behind the front is uniquely defined by the 
state ai in front of the front [see Fig. f|(b)]. Hence there 
is functional relation between them: 



a/ = A(a,i). 



(27) 



Together with ([23|) we have a uniquely defined relation 
between c, a^, and a/. This behavior fits into the general 
picture of front propagation in bistable systems |l9[ (see 



also Sec. [V B| ). If we know the function A we can cal- 
culate the values of ax and a 2 by solving A -1 (a) = A(a) 
[see also Fig. |](b)]. Note that the values of ai and a 2 are 
in general irrational [^6| . The value of a determines only 
the sizes of the domains. Whether the chain is commen- 
surate or incommensurate is irrelevant. But numerically 
we have never found such states for values of a near in- 
teger multiples of 2ir. 
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FIG. 6. Velocity-force characteristics for the uniform slid- 
ing states and the two-domain states. Thick (thin) lines 
denote two-domain (uniform sliding) states. Solid (dotted) 
lines indicate stable (unstable) solutions. The parameters are 
N = 2584, M = 987, b = 2, and 7 = 0.5. 

The selected values of a\ and a 2 are of course functions 
of the applied force F. Figure || shows the velocity- force 
characteristic of the two-domain states. By varying F 
they disappear because of two reasons. First, one of the 
domains shrinks to zero and the velocity-force character- 
istics of the two-domain state and uniform sliding state 
come together. Secondly, a\ or a 2 may move onto an 
unstable branch of the velocity-distance characteristic. 
Thus, the two-domain state still exists but it has become 
unstable. In Fig. || the velocity-force characteristic of 
this unstable two-domain state is denoted schematically 
by dotted extensions. 



B. Quasi-continuum description of the front 

We have seen that the state of the domains can be 
described by (7|). They are characterized by a\^ 2 and 
V\,i- To describe the fronts in the same way we assume 
that a and v are continuous functions which are varying 
slowly in space and time. The space coordinate is the 
particle index j. It becomes a real variable in a quasi- 
continuum description. Thus, we write Xj(t) = x(j,t). 
The discrete Laplacian + 2^-1 — can be written 
as an infinite series of differential operators, i.e., 



x j+1 



(t) + xj-i(t) - 2xj(t) = 4sinh 2 Qc^ x(j,t). (28) 



We generalize the ansatz (0) by assuming that <p is a 
function of j and t, i.e., 
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5CJ.*) = V0*.0 + /(V0'.*)), f{<P + 2tt) = (29) 



where / is the solution of the hull function equation 
Now we define local values of a and v by 



<9j<y9, and 



(30) 



Plugging the ansatz ( p9| ) into the equation of motion (Q) 
and averaging over the phase (p we get 



5 4 w = J D(a j )9 7 a + F-F a (a,«), 
dta = djV, 



(31a) 
(31b) 



where 



D{x) 



( sinhx/2 
V x/2 



l + -+0(x% (31c) 



and Fjj(a, v) is the velocity-force char acte ristic of the uni- 
form sliding state given by ([H]). Eq. (31a) is only approx- 
imately correct because we have assumed that the hull 
function / does not depend on a and v. Furthermore, 
the ansatz (^9|) cannot be exact in a front. Nevertheless, 
the approximation (|3l]) is correct in leading order of a 
multiple-scale perturbation theory |2^J28|]. 




FIG. 7. Schematical drawing of the velocity-distance char- 
acteristic and the corresponding nonlinearity g. 

In general eq. ( |3l| ) cannot be solved analytically. But 
we are able to discuss the front solutions qualitatively. 
Assuming stationarity of the front in the co-moving frame 
we get 



a(j, t) = a ti ~ cb), v(j, t) = v{j - ct), 



(32) 



where c is the front velocity (|23]). From ( 31b ) we get 
—ca' — v 1 which can be integrated leading to 



Plugging fl3J) and (g3j) into 
first two terms of D yields 

1 



la|) and keeping only the 



(1 - c )a' + —a!" + g(a, v , c) = 0, 



with 
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g(a,v ,c) = F - Fu(a,v - ca). 



(34) 



(35) 



Eq. (|33|) means a straight line in the velocity-distance 
characteristics (see Fig. Front solutions exist for those 
values of vq and c for which (|3^) intersects the velocity- 
distance characteristic of the uniform sliding state three 
times. The two outside intersection points have to lead 
to stable uniform sliding states (characterized by a\ and 
CL2) whereas the inner point has to belong to an unstable 
uniform sliding state. This is the reason why two-domain 
and multi-domain solutions appear only near resonance 
points where the velocity-force characteristic has a nega- 
tive slope (see Fig. ^|). Because there is no resonance for 
a/27r integer, we understand why we have not found two- 
domain solutions for values of a close to integer multiples 
of 2tt. 

If the requirements on vo and c are fulfilled, the non- 
linear term g in (p4J) will have three nodes and will be 
N- shaped (see Fig7]7|). A front solution is a heteroclinic 
orbit of (BjO) which goes from a\ to or vice versa. Thus, 
we are looking for solutions with boundary conditions 
o(— 00) = ai,2 and a(oo) = ai,\- A heteroclinic orbit oc- 
curs only if the unstable manifold of the fixed point a\ t 2 
is the stable manifold of 0,2.1- This is possible only on 
a one-dimensional manifold in the parameter space of vq 
and c. Thus ( gjj ) is justified. 

To calculate the stable and unstable manifolds we lin- 
earize (|3~i| ) around the fixed points 01,2- For the pertur- 
bation 5a = a — a± 7 2 we make the ansatz 6a = cxp(Aj) 
which leads to the characteristic polynomial 



— A 3 + (l-c 2 )A + a a5 ( ai , 2 ) =0. 



(36) 



(33) 



Because of <9 a fl(ai.2) > (see Fig. ^) there is one nega- 
tive root Ai < 0. If (1 - c 2 ) 3 + {3d a g/A) 2 > the two 
other solutions are conjugated complex with a real part 
that is just — Ai/2. Numerically we always found sub- 
sonic front motion (i.e., \c\ < 1) leading to an unstable 
manifold that spirals out of the fixed point. Thus, the 
precursor of the front is non-oscillatory whereas its tail is 
oscillatory because the particles have inertia and there- 
fore respond with an exponentially decreasing oscillation 
after an acceleration or deacceleration. 

In order to verify this qualitative picture numerically 
one has to extract a and v from the data. In principle this 
could be done by local fits of the dynamic hull function 
from which we obtain tp(j, t) and subsequently a and v. 
But this is a very tedious way which is not necessary be- 
cause we are only interested in the qualitative form of the 
shape of the front. The following simple method is suffi- 
cient for this task. For uniform sliding states it leads to 
values of a and v which are identical with the exact ones. 
We introduce for each particle a sequence of times t n j 
defined by Xj(t n j) = (2n + l)ir and ij(t n _j) > 0. From 
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these sequences we get the following approximations of v 
and a: 



2n 



"7l,j vfl — ljj 



-, a(j,t) & v(j,t)(t 



m,j — l ^71, J ) •> 



(37) 



where n and m are chosen in such a way that (i) t n -i,j < 
t < t n _j and (ii) t TO)J _i is the time closest to i.e., 
|£mj-i — t n j\ = mina \taj-i — t n .j\. This definition of a 
is necessary in order to avoid spurious values which differ 
from the expected value by ±27r. 



250 
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FIG. 8. Quasi-continuum description of the two-domain 
solution. The values of v and a are obtained from the formulas 
(fi7|). 20 snapshots taken at different times and shifted by ct 
are superposed in order to get the details of the fronts. The 
value of c is chosen in such a way that the superposition yields 
a smooth curve. The best value of c is 0.516. The parameters 
are the same as in Fig ^. 

For a fixed value of t one gets a snapshot of v and a. 
The superposition of many snapshots shifted by ct gives 
the impression of a smooth curve (see Fig. ||). We have 
tuned c until the curves are as smooth as possible. It 
turned out that this method is a very accurate way to 
measure the front velocity c. The precursors and tails 
of both fronts are non-oscillating and oscillating, respec- 
tively. Furthermore, the oscillatory tail of the front on 
the left decays roughly two times slower than the precur- 
sor of the front on the right. Both observations are in 
full agreement with our analytical reasoning above. 



C. Multi-domain states 

Starting from an arbitrary initial condition one gets 
either a uniform sliding state (if a stable one exists) or a 
multi-domain state but only rarely a two-domain state. 
This is especially true for large systems. It is a very gen- 
eral behavior of bistable spatially extended systems, at 
least for the initial phase of the dynamics. Different parts 
of the system establish themselves independently into one 
of the bistable states. It is therefore natural that several 
domains occur. After the initial formation of a multi- 
domain state, domains may shrink and eventually disap- 
pear on a slow time scale. This can be understood by the 
fact that an attractive force between the fronts exist [^9) . 
This force is caused by the overlap of the precursor and 
the tail of the neighboring fronts. It can be calculated by 
singular perturbation theory |l9|]29| ]. In the case of non- 
oscillating precursors and tails the force decreases expo- 
nentially with distance (^9|. In our case where the tail is 
oscillating, the resulting force is also oscillating |]l9| , ^o| . 
Therefore equilibrium positions are possible where a pair 
of two fronts form a bound molecule. The available dis- 
tances between the fronts are quantized. Figure [9] shows 
an example where this quantization is clearly seen. In 
accordance with Shilnikov's theorem such molecules are 
possible only if the oscillatory tail decays slower than the 
non oscilla tory one [ plp^ ] which is indeed the case (see 
Sec. IVB ). As a consequence the system shows spatial 
chaos. 




100 150 200 

j-ct 
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FIG. 9. An example of a multi-domain state. The param- 
eters are the same as in Fig ^. 

For values of F near higher order resonances the 
velocity-distance characteristic of the uniform sliding 
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state shows several resonance folds. Thus, multi-domain 
states are possible which are built up from more than two 
different domain types. All our numerical experiments 
have shown that no matter how many domain types oc- 
cur in the transient, at the end (i.e., in the long-time 
limit) only two or three domain types survive. Further- 
more, all fronts travel with the same velocity c. There 
is a simple argument why more than three domain types 
are inconsistent with the last fact. Consider the case of 
four different domain types with <zi < 02 < 03 < 04 and 
V\ > i>2 > U3 > V4. Between two consecutive values <2; 
and di+i there should be no additional stable state and 
only one unstable state. Thus, if two of such domain 
types are neighbors the states are functionally related in 
accordance with (p7|). Let us assume a sequence of do- 
mains from rig ht (domain type 1) to left (domain type 
4). Using ( P7| ) we get a i+ i = A{ai). This sequence is 
therefore uniquely determined by ai . Now the condition 
that all fronts have the same speed cannot be fulfilled 
because there is only one free variable but (at least) two 
conditions, namely, C12 = c 2 3 = c 34 . Thus, four different 
domain types are in contradiction with the observation 
concerning the front velocities. 

For F < b domain-like solutions occur where in one 
domain state the particles sit in potential wells (aver- 
age particle distance a is an integer multiple of 2tt, e.g., 
a = 0) and do not move. Such so-called traffic-jam solu- 
tions have been already reported by Braun et al. jlB] |l7) . 
They mainly occur near integer values of a/2-K where the 
domain-like states discussed in this section are not pos- 
sible. 



V. CONCLUSION 

In this paper the regular sliding states (periodic and 
quasi-periodic solutions) of the FK model are investi- 
gated for large chains (N > 100) and for values of 7, 
&, and F where all forces in the equation of motion (|l|) 
are of the same order. 

Instead of classifying these attractors to be either pe- 
riodic or quasiperiodic, a more informative distinction is 
whether their locally average particle density is uniform 
or not. In a periodic state the density is uniform be- 
cause all particles perform the same motion only shifted 
in time. The motion of the whole chain is completely 
determined by a single periodic function, the dynamical 
hull function. The velocity-force characteristics of the 
periodic attractors show peaks at certain values of the 
velocity. These peaks are caused by resonances of the 
"washboard" wave. In the frame of the center of mass 
of the chain the external potential can be seen as a wave 
(washboard wave) with wave number a (which is the in- 
verse average particle density) and frequency v (which is 
the average sliding velocity) . Resonances occur for those 
values of v where n washboard waves are able to decay 
into a phonon [wavenumber k, frequency u)(k)] in accor- 
dance with "momentum" and "energy" conservation [i.e., 
k = na and u>(k) = nv]. A similar process with a decay 
into two phonons (parametric resonance) explains why 
the uniform sliding state may be unstable even though 



the differential mobility is positive. 

Quasiperiodic states emerge from instabilities of uni- 
form sliding states. They are characterized by domains of 
different average particle densities. In each domain the 
average particle velocity is uniform and nonzero. The 
walls between neighboring domains move all with the 
same velocity. This domain-wall or front velocity is dif- 
ferent from the velocity of the center of mass. Thus, 
each particle goes through all domains. In a quasiperi- 
odic attractor not more than three different domain types 
are possible. For chains close to the most commensurate 
case (i.e., a is an integer multiple of 2n) the average par- 
ticle velocity is zero in one domain. In this so-called 
traffic-jam state |H| the particles are alternately switch 
between stationarity and sliding. Because the traffic-jam 
state exists also for finite temperatures |15 ItJ states with 
different sliding domains will presumably survive finite 
temperatures. 

In this paper we have developed a quasi-continuum de- 
scription of the FK model on the basis of slowly varying 
locally averaged inverse density a and velocity v. That 
is, a and v slowly depend on the particle index j and the 
time t. A set of partial differential equations [Eqs. J3l|)] 
governs the dynamics of the variables (assuming j to be 
a continuous space variable). It is not possible to solve 
these partial differential equations analytically. But it 
is quite helpful to understand (i) why the FK model or- 
ganizes itself into domain-like states with domain states 
which are not determined from outside, (ii) why there 
are not more than three different domain types, (iii) why 
the size of the domains is quantized, (iv) why states with 
several domains of the same type are possible, (v) why 
a multi-domain state can be seen as an example of spa- 
tial chaos, and (vi) why for a fixed value of the external 
force F the number of stable states increases exponen- 
tially with N. The huge number of stable states leads 
to multi-hysteretic behavior like in a ferromagnet where 
the position in the velocity-force characteristic strongly 
depends on the history. 

The phenomenon of domain-like states where each do- 
main is characterized by a spatially periodic but station- 
ary solution has already been found in hydrodynami- 
cal pattern formation (Taylor-Couette system |33| and 
Rayleigh-Benard system |Q). Theoretically this be- 
havior has been modeled by nonlinear phase equations 
@|5|j36|. In our case the domain states are periodic in 
space and time, they are traveling waves. 

In experiments where locally resolved measurements 
are not possible, the most important consequence of the 
domain-like sliding states are (i) quasiperiodicity in the 
time signals like the velocity of the center of mass, (ii) 
flattening of the resonance peaks (see Fig. ^|), and (iii) 
multihysteretic behavior. Many of these features should 
disappear for small values of N if they are caused by 
domain-like states. Because the inertia term is important 
for these states, we do not expect it in CDW-systems. 
In adsorbate layers and ionic conductors the appearance 
of them should be possible but it may be difficult to 
drive them uniformly and to measure the velocity-force 
characteristic. The ideal system to check our theory are 
Josephson-j unction arrays because the force and the ve- 
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locity correspond to the driving current and the voltage, 
respectively. Furthermore, the damping constant 7 and 
the number of junctions N can be chosen by the fabrica- 
tion process. The average distance a — 2irM/N is also 
easily accessible because the number M of flux quanta 
can be chosen by the initial preparation of the system. 
Thus, for rings of more than 50 Josephson junctions we 
predict the occurrence of domain-like sliding states. 
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powers of two. Now we are able to calculate {Ii, . . . , Im} 
from {fx, ... , /m} by two Fast Fourier Transformations 
(FFT). The first one is an inverse FFT which calcu- 
lates the hull function in real space. Next we calculate 
sin[v5 + /(<p)]. The second FFT yields {h, I M }. 

In order to get reliable results the value of the cutoff 
M has to be chosen carefully. Due to the FFT the hull 
function is approximated on a lattice with lattice con- 
stant n/M. Because of (0) this corresponds to a time 
resolution At = ir/(vM) oixj(t). The fastest time scale 
in the system is given by 2ir divided by the maximum 
phonon frequency which is 2 in accordance with fll4|). 
Assuming that the fastest time scale has to be resolved 
at least by two steps At we get reliable results only if 



> 



M' 



(A5) 



APPENDIX A: NUMERICAL AND 
ANALYTICAL APPROXIMATIONS OF THE 
HULL-FUNCTION EQUATION (|) 

In this appendix we describe the numerical scheme we 
have used to solve the differential delay equation (^|) for 
the dynamic hull function f(<p). In the simplest case this 
scheme leads also to a nonlinear analytical approxima- 
tion. 

Because of (0) and (|o|) we expand the dynamic hull 
function / into a Fourier series 



00 

E 



fn 



im<p 



+ C.C. 



(Al) 



In the numerical approximation we replace this infinite 
series by a finite one with a cutoff M, i.e., J2m=i ~ * 
Sm=i- P m g§ m S this ansatz into the differential delay 
equation (^) leads to a set of nonlinear algebraic equa- 
tions for the coefficients {/1, / 2 , ...,/,„,.. .}: 

|(™) 2 - w 2 (ma) - ijmv]f m = W m (/i, /2, • • • , fj, ■ ■ •)> 

(A2) 

where w(k) is the phonon dispersion relation ( |l^ ) and 



Im(fl,fc, ■■ ■ ,fj, ■■ •) 



2vr Jo 



(A3) 



The driving force F does not appear in the set of equa- 
tions ([A2]). That is, we first obtain a solution of (A2) for 
a given value of v. After that we get the corresponding 
value of F by using (|l2|), i.e., 



All results reported in this paper are obtained for M = 
32. Thus, v has to be larger than 0.06. 

For M = 1 we get an analytic result parameterized by 
the modulus of /1 . We write 



fi 



A 



(A6) 



Due to the integral representation of the Bessel functions 
of the first kind we can obtain 1\ (/1). Together with ( [A2] ) 
we get 

[v 2 - w 2 (a)]A = b[J 2 (A) - J (A)} sin^, (A7) 
jvA = b[J 2 (A) + J (A)] cos tp. (A8) 

The elimination of ip leads to polynomial of second order 
in v 2 . Thus, we get an analytic solution parameterized 
by A: 



v(A) = 



\ 



7 2 



\ 



with 



and 



Jo(A) - J 2 {A) 
J a (A) + J 2 (A) 



b=[J (A)-J 2 (A)]b. 



Using (A4) we get 



F(.ll =-,■(.!) (1 + ^- 



(A9a) 



(A9b) 



(A9c) 



(A9d) 



F = 1 v{\ + 2y j m 2 \f m \ 2 ) 



(A4) 



In the numerical approximation we solve the set of M 
algebraic equations ( |A2| ) with the Newton method. The 
most time-consuming part of the computation is the cal- 
culation of I m . In order to speed up the computation 
we restrict ourselves to values of the cutoff M which are 



APPENDIX B: LINEAR STABILITY ANALYSIS 
OF THE UNIFORM SLIDING STATE 

In this appendix we explain our procedure to calcu- 
late numerically the stability of the uniform sliding state. 
Furthermore, we derive an approximation for v c . 



10 



Plugging the Floquet-Bloch ansatz ( p^ ) into the lin- 
earized equation of motion (|l^) leads to 

v 2 4&) + (2A fe + l)vc' k {ip) + (X 2 + A fc7 )cfeM = 

c k (tp ~ a)e^^r- k + c k (ip + a)e^ fe - 2c k {p) 
-bcc>8[<p + f(ip)]ck(<p). (Bl) 

The Fourier ansatz 



where 



(<p) = J2 c k , m e mifi + c.c. 



(B2) 



m=0 



turns this equation into a set of infinitely many linear 
equations: 



Li." j ma + ) + (Afc + imv) 2 + ^y(\ k + imv) 



Ck- 



-b ]T ( K ™ 

—m'Ck.m' 

+ K m 

+m' c k,m' ), (B3) 



where 



1 f 2 " 

K m = — / cos[ip + f(<p)]e 
27r Jo 



(B4) 



Again we have to choose a cutoff M' to solve this set 
numerically. In order to be consistent with the cutoff 
M = 32 of the Fourier expansion of the hull function / 
we have chosen M' = 15. Because of ( |l9|) the stability 
depends on the number N of particles. For large N the 
eigenvalues A^ lead to a continuous function X(2Ttk/N). 
Since we are mainly interested in large systems (N > 50) 
we have chosen N — 100. It is difficult to check numer- 
ically whether the uniform sliding state is stable or not 
because of the Goldstone mode Ao = 0. Our numeri- 
cally obtained value of Ao fluctuates around zero because 
of unavoidable errors. Therefore our instability criterion 
reads: The uniform sliding state is unstable if at least one 
eigenvalue is larger than 0.05. This value is considerably 
larger than the amplitude of the numerical fluctuations 
of A . 

We are able to obtain an analytical approximation of 
the largest sliding velocity v c at which parametric reso- 
nance is just able to destabilize the uniform sliding state. 
We do this for zero damping because physical intuition 
tells us that v c decreases monotonically with 7. In fact 
for 7 <C 1 and b = 0(7°), the largest sliding velocity is 
the undamped one minus a correction term of order 7. 
The approximation makes three assumption: (i) Para- 
metric resonance at v c occurs for that value of q which 
maximizes vf(q). For < a < tt this implies q = tt. 
(ii) All Fourier coefficients of Ck are zero except Ck,o and 
Ck,\- (iii) The integrals K m are calculated only in leading 
order of b which yields 



K„ 



2tt 



Sl.m + 6-1,: 



(B5) 



where 5 n _ rn is the Kro necker symbol. With 7 = and 
these assumptions (B3) reduces to 



, 2 + A 2 
6/2 



6/2 

(\ + iv) 2 



Cfc,0 
Cfc, 1 



(B6) 



lo = uj(k) = u>(a + k) = u)(a/2 ± tt) = 2 cos . (B7) 



A nontrivial solution of (B6) implies a zero determinant 
leading to a characteristic polynomial of second order in 
(A + it>/2) 2 . Solutions with nonzero real part of A occur 
only if (uj 2 - v 2 /A) 2 < 6 2 /4. Therefore we get (2l]). 
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